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Abstract 

In this paper, we consider multiple signals sharing same instantaneous frequencies. 
This kind of data is very common in scientific and engineering problems. To take ad¬ 
vantage of this special structure, we modify our data-driven time-frequency analysis by 
updating the instantaneous frequencies simultaneously. Moreover, based on the simul¬ 
taneously sparsity approximation and fast Fourier transform, some efficient algorithms 
is developed. Since the information of multiple signals is used, this method is very ro¬ 
bust to the perturbation of noise. And it is applicable to the general nonperiodic signals 
even with missing samples or outliers. Several synthetic and real signals are used to 
test this method. The performances of this method are very promising. 


1 Introduction 

Data is an important bridge connecting the human being and the natural world. Many 
important information of the real world is encoded in the data. In many applications, fre¬ 
quencies of the signal are usually very useful to reveal the underlying physical mechanism. 
Therefore, in the past several decades, many researchers devoted to find efficient time fre¬ 
quency analysis methods and many powerful methods have been developed, including the 
windowed Fourier transform m , the wavelet transform mm, the Wigner-Ville distribu¬ 
tion [9], etc. Recently, an adaptive time frequency analysis method, the Empirical Mode 
Decomposition (EMD) method 115], [20] was developed which provides an efficient adap¬ 
tive method to extract frequency information from complicate multiscale data. Inspired by 
EMD, recently several new time frequency analysis methods have been proposed see e.g. the 
synchrosqueezed wavelet transform [6], the Empirical wavelet transform [B], the variational 
mode decomposition [10]. 

In the last few years, inspired by the EMD method and compressive sensing 131 El E], 
Hou and Shi proposed a data-driven time-frequency analysis method based on the sparsest 
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time-frequency representation [12]. In this method, the signal is assumed to fit following 
model: 

K 

f(t ) = y^a k (t) cos 6 k (t), t 6 [0,1], (1) 

fc=i 

where ak(t), 9k(t) are smooth functions, 0' k (t ) >0, k = 1, • • • ,M. We assume that ak(t) 
and 6' k are ’’less oscillatory” than cos 9kit)- After some scaling, we can always make the 
signal lies in the time interval [0,1]. So, here we assume the time span of the interval 
is [0,1]. We borrow the terminology in EMD method and also call ak(t) cos 9k(t) as the 
Intrinsic Mode Functions (IMFs) [T5j. This model is also known as Adaptive Harmonic 
model (AHM) and is widely used in the time-frequency analysis literatures [1 [ [6]. 

The AHM model seems to be natural and simple In this model, the main difficulty is 
how to compute the decomposition. In the AHM model, the number of degrees of freedom 
is much larger than the given signal. This makes the possible decomposition is not unique. 
One essential problem is how to set up a criterion to pick up the ’’best” decomposition 
and this criterion should be practical in computation. Inspired by the compressive sensing, 
we proposed to decompose the signal by looking for the sparsest decomposition. And the 
sparsest decomposition is obtained by solving a nonlinear optimization problem formulated 
as following: 

K 

min K, Subject to: / = Y"' cos 9k, ak cos 9k &T>. (2) 

( fl fc)i<fc<jf ,(h)i<fc<M k _^ 

where T> is the dictionary consist of all IMFs which makes the decomposition satisfies the 
AHM model (see [T2] for precise definition of the dictionary). 

Two kinds of algorithms were proposed to solve the optimization problem ©• The first 
one is based on matching pursuit H3H2] and the other one is based basis pursuit mm- 
The convergence of the algorithm based on matching pursuit was also analyzed under the 
assumption of certain scale separation property [141 . 

In the previous works, we consider the case that only one signal is given. However, 
in many application, we can obtain many different signals and these signals have same 
frequency structure. For instance, in the monitoring of buildings, usually many sensor are 
put in different locations of the same building to measure the vibration. From these sensors, 
many measurements of the vibration are obtained. Since the signals from different sensors 
measure the same building, they should have same instantaneous frequencies which associate 
with the natrual frequencies of the building. If we analyze the signals from different sensors 
individually, the structure that these signals share the same instantaneous frequencies is 
wasted. Intuitively, by taking advantage of this structure, we may get more robust and 
efficient methods. In this paper, we will propose several methods by exploiting this special 
frequency structure. 
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In this paper, we consider multiple signals sharing same instantaneous frequencies. Nat¬ 
urally, we modify the adaptive harmonic model ([!]) to deal with this kind of signals. 

K 

f{t) = ^2 cosd k (t), t € [0,1], J = 1, • • • , M, (3) 

k =1 

where M is the number of signals, f J is the jth signal. For each j = l,--- , M fixed, it 
is the AHM model 0, i.e. 4 (t), &k(t) are smooth functions, 0' k (t) >0, k = 1, ■ ■ • ,K. 

We assume that a k (t) and 9' k are ’’less oscillatory” than cos 9k(t). The main feature of 
this model is that different signals have same phase functions 9k- We call (0 Multiple 
Adaptive Harmonic model (MAHM). 

Inspired by our previous work on single signal, we also look for the sparsest decompo¬ 
sition which satisfy the MAHM model 0 by solving following optimization problem: 

min K (4) 

( n j a 

K a k"k)k=t,... ,K 

K 

Subject to: f\t) = cos 9k(t), a J k (t) cos 9k(t) 6 V, k = 1, • • • , K, j = 1, • • • , M. 

k =l 

The dictionary V we used in this paper is same as that in [T2j • It can be written as 

V = {acos9: a g V(9, A), 9' € V(9, A), and 9\t) > 0, Vf € [0,1]} . (5) 

V(9. A) is the collection of all the functions ’’less oscillatory” than cos 9(t). In general, it is 
most effective to construct V(9, A) as an overcomplete Fourier basis given below: 

V(9,X) = spanjlYcos 1’ 

[ V \ ZL e J J i<k< 2 \Lg V \ ZL eJJ i<k< 2 XLg) 

where Lg = L'J i s the largest integer less than (•), and A < 1/2 is a parameter 

to control the smoothness of V(9, A). In our computations, we typically choose A = 1/2. 

In the rest of the paper, we will introduce several algorithms to approximately solve 
above optimization problem 0. First, we give a generic algorithm based on matching 
pursuit and nonlinear least squares in Section 2. This algorithm can be accelerated by 
fast Fourier transform (FFT) if the signal is periodic. This will be reported in Section 3. 

For general nonperiodic signals, in Section 4, we also develop an efficient algorithm based 
on group sparsity and the algorithm in Section 2. This algorithm can be also generalized 
to deal with the signals with outliers or missing samples (Section 5). In Section 6, we 
present several numerical results including both synthetic and real signal to demonstrate 
the performance of our method. At the end, some conclusion remarks are made in Section 
7. 
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2 Method based on matching pursuit 

Following the same idea as that in [12j . we proposed a method based on mathcing pursuit 
to get the sparsest decomposition, see Algorithm [I] 

Algorithm 1 Method based on matching pursuit 
Input: Signals j = 1, • • ■ , M. 

Output: Phase functions and the corresponding envelopes: 9 k , a k , k = 1, ■ ■ ■ , K, j = 

1 M. 

1: Set r{ = j = 1, • • • , M and k = 1. 

2: while max ^||r^||; 2 ^ > eo do 

3: Solve the following nonlinear least-square problem: 

M 

min ^ || r J k - a{ cos 6 k \\p (7) 

a k$k j = l 

Subject to: 0' k > 0, a k (t) € V(9k), j = 1, • • • , M,. 

4: Update the residual 

r J k+ 1 = r[ - a{(t) cos9 k , j = 1, • • • , M. 

5: Set k = k + 1. 

6: end while 


To solve the nonlinear least-squares problem ([T]) in Algorithm [TJ we use the well-known 
Gauss-Newton type iteration. This algorithm is also very similar as that in [T2j. The only 
difference is the update of the phase function. In our model, different signals share the same 
phase function. So in the algorithm, we only update one common phase function by taking 
some average among different signals. 

By integrating Algorithm [l] and [2] together, we can get a complete algorithm to compute 
the sparsest decomposition. The key step and also most expensive step is to solve the least- 
squares problem ([8]). It does not need too much time to solve ([8]) once, however, we have to 
solve this least-squares problem many times to get the final decomposition. This makes that 
the algorithm is not very practical in applications. In next two sections, we will propose 
some acceleration algorithms. 
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Algorithm 2 Gauss-Newton iteration to solve the nonlinear least-squres problem 
Input: Initial guess of the phase function 0 k = Oq. 

Output: Phase functions and the corresponding envelopes: Ok, a J k , j = 1, • • • , M. 
1: while || 0 k +l - 0 'l\\i 2 > e 0 do 

2: Solve the following least-square problem for each signal r{_ 1 , j = 1, • • • , M: 

^n+^ln+i H “ °fc n+1 W cos0 fe(O “ ? 4’ n+1 ( i ) sin 0^(4) ||p 
a i’ 

Subject to b{’ n+1 (t) G V{0 k ). 


3: Calculate averaged updation of frequency: 


A Wj = 


>f +1 (r +I )'-<4'" +1 (»r +1 )' 

-^-A-V -A Au = 

( 4 ” +1 ) +(* 4 ’” +1 ) 


^ = 1 _ (Q) 


Z^j=l 1 k 


where r£ n+1 = ^4’ n+l ) + ^’” +1 j and (•)' denote the derivative of (•) with respect 


4: Update 0 k 


A 6 = [ Aco(s)ds, 0l +1 = 01- PAO, 
Jo 


where f3 G [0,1] is chosen to make sure that 0 k +l is monotonically increasing: 


/3 = max < a G [0,1] : — (0 k — aAO) > 0 


5: end while 
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3 Periodic signals 


If the signals are periodic, we can use a standard Fourier basis to construct the V(6, A) 
space instead of the over-complete Fourier basis given in (Jb]). 

V P (6,X) = spanilYcos (j^X) ’( sin (j~)] ( 10 ) 

[ V \ L eJ J i<k<XL e V \ L dJ J \<k<\L e J 

where A < 1/2 is a parameter to control the smoothness of functions in V p (9, A) and Lq = 
(0(1) — 0(O))/27t is a positive integer. 

In this case, we developed an efficient algorithm based fast Fourier transform to solve 
the least-squares problem Q in p[2j. To make this paper self-contain, we also include the 
algorithm here. 

Suppose that the signal r k - 1 is measured over a uniform grid tj = j/N, j = 0, • • • , N— 1. 
Here, we assume that the sample points are fine enough such that r k ~i can be interpolated 
to any grid with small error. Let 6 = g^Jg^Q) be the normalized phase function and 

Lg = which is an integer. 

The FFT-based algorithm to approximately solve ([8]) is given below: 

Step 1: Interpolate rk- 1 from {U}^L 1 in the physical space to a uniform mesh in the (In¬ 
coordinate to get rgn and compute the Fourier transform rg«: 

r 0£, j = Interpolate (r k -i,9%j) , (11) 


where 9 k j' J = 0, • • • , IV — 1 are uniformly distributed in the ^-coordinate,i.e. j = 

2ttLq j/N. And the Fourier transform of rg« is given as follows 

1 N 

= “ = —W/2 + 1, • • ■ ,N/2, (12) 

3 = 1 


, TV S k i~ e k 0 

where 9 kJ = & L • . 


Step 2: Apply a cutoff function to the Fourier Transform of rgn to compute a and b on the 


mesh of the ^/-coordinate, denoted by ag« and bgn: 

k k 


a 9 ]] 


= T 


—l 


—l 


rgn (uj + Tgr^ + fgn (u - Lgn'j'j ■ X \ (u/Lffn'j 

rei (w + Lgn ) - rgn (u ; - Tg» ) ) • xx (w/Tg? 


bgn = J- 


T 1 is the inverse Fourier transform defined in the 6 coordinate: 


T 


—l 


re " ' = N 


N/2 

£ 

u=-N/ 2+1 


rgne 


i2-KU]8j. 


, j = 0, ■ ■ ■ , N - 1. 


(13) 

(14) 


(15) 
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Step 3: Interpolate ag and bg from the uniform mesh { 9 £ j}jLi in the ^-coordinate back to 
the physical grid points 

a(U) = Interpolate (ag,U) , z = 0, ■ ■ ■ , IV — 1, (16) 

b(U) = Interpolate ( bg,U ), i = 0, • ■ ■ , N — 1,. (17) 


The low-pass filter xa(w) in the second step is determined by the choice of V p (9, A). 
In this paper, we choose the following low-pass filter xa(w) to define V p (9, A): 


XaM 


1 + cos(7to;/A), —A < oj < A 
0, otherwise. 


(18) 


By incorporating the FFT-based solver in Algorithm [2j we get an FFT-based iterative 
algorithm which is summarized in Algorithm [3} 


Remark 3.1. We remark that the projection operator P\/ v (0-, g) /act a low-pass filter in 
the 9-space. For non-periodic data, we apply a mirror extension to A 9 before we apply the 
low-pass filter. 

Gauss-Newton type iteration is sensitive to the initial guess. In order to abate the depen¬ 
dence on the initial guess, we gradually increase the value ofr] to improve the approximation 
to the phase function so that it converges to the correct value. The detailed explanation can 
be found in TWf . 


4 nonperiodic signals 

In most of the applications, the signals are not periodic. To apply the FFT-based algorithm 
in previous section, we have to do periodic extension for general nonperiodic signals. In 
this paper, the signals are assumed to satisfy the MAHM model ([3]). One consequence of 
the MAHM model is that the phase function can be used as a coordinate and the signals 
are sparse over the Fourier basis in 9 coordinate. Then, one natural way to do periodic 
extension is to look for the sparsest representation of the signal over the over-complete 
Fourier basis defined in 9 coordinate. 

More precisely, the extension is obtained by solving an l\ minimization problem. 

min 11 x 11 1 , subject to <&g • x = f (20) 

xeM^ 

f is the sample of the signal at tj,j = 1, • • • ,N S , N s is the number of sample points, tj 
may not be uniformly distributed, however, we assume that the sample points are fine 
enough such that f can be interpolated over any other grids without loss of accuracy. 

G C N ° xNb is a matrix consisting of basis functions, Ng is the number of Fourier modes. 
*o(j,k) = e ikO(tj)/L e , j _ 1 ; ... ; JV S , k = —Ng/2 + 1, ■ ■ ■ , Nb/2. Lg is a positive integer 
determined by the length of the period of the extended signal which will be discussed later. 
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Algorithm 3 (FFT-based algorithm to solve the nonlinear least-squres problem) 

Input: Initial guess of the phase function Of = Oq, 7] = 0. 

Output: Phase functions and the corresponding envelopes: Ok, aj., j = 1, • • • , M. 

1: while rj < A do 

2: while \\0f +1 - 0f\\ P > do 

3: Interpolate rk-\ to a uniform mesh in the ^-coordinate to get rgn and compute 

the Fourier transform rgn. 

k 

4: Apply a cutoff function to the Fourier Transform of rgn to compute a and b on the 

mesh of the ^-coordinate, denoted by agn and bgn. 

k k 

5: Interpolate ag and bgn back to the uniform mesh of t. 

6: Calculate averaged updation of frequency: 


Aujj = 


4" +1 (r +1 )'-r +1 (°r +1 )' 

(< + 1 ) 2 +( h ”* +1 ) 2 


, At o = 


LfU rf +1 

2^7=1 1 k 


where r ),’ n+1 = ^a^ n+l ^ 1 'j and (•)' denote the derivative of (•) with 

respect to t. 

Update Of, 

AO’ = Pv(e-, V ) (Aw), A 0(t) = f A 0\s)ds, 0f +l = Of - /3A9, 

Jo 

where f3 € [0,1] is chosen to make sure that 0 f +1 is monotonically increasing: 


(i = max | a G [0,1] : — (Of — a AO) > 0 . 

and P Vp{ 0 . v) is the projection operator to the space V p (0: rj). 

end while 

i] = ij + Ar] 

end while 





To get the periodic extension of M signals, one way is to solve (20) M times inde¬ 
pendently. However, these M signals are not independent. They share the same phase 
functions. To take advantage of this structure, inspired by the methods for simultaneous 
sparsity approximation nansi, we propose to solve following optimization problem to get 
the periodic extension of M signals simultaneously. 


min IIXI 

xeK jV 6 xM 


2 , 1 ) 


subject to • X = F 


( 21 ) 


where X is a N b X M matrix, X = {xj^)j=i,---,N b) k=i,-~Mi F = [f 1 , - - - , f M ] is a N s x M 
matrix and 

N b / M \ 1 / 2 


|X|| 2 ,i = 

j =i Vfc=i 


( 22 ) 


The problem left is how to solve the optimization problem (21) efficiently. Since the sample 
points are fine enough, we can assume that the sample points are uniformly distributed in 
6 coordinate, i.e., 0(tj) = (j — 1)A 0 and A 6 = 2TrLg/Ng. Otherwise, we could use some 
interpolation method (for instance cubic spline interpolation) to get the signals over the 
uniformly distributed sample points. Now, we define a matrix, Qq E C NbXNb , consisting of 
complete Fourier basis, 


®e{j,k) = 


1 


Jk(j-1)A6/L e 


VN~b 


j = !,■•• ,N b , k = !,••• ,N b . 


(23) 


&g(j, k) represents the element of at jth row and kth column. Using the property of 
the Fourier basis, we know that is an orthonormal matrix which will give us lots of 
convenience in deriving the fast algorithm. 

For the convenience of notation, denote H be the set of index of F, 


n = {(j,k):j = !,-■■ ,N b , (j - 1)A0 < 0(1) - 0 ( 0 ), k = ,N b }. 


(24) 


Then, for any (j. k) E H, F(j, k) is a sample point. We also define an extension of F by zero 
padding, 

F (j,k), (j,k)en, 

0, otherwise. 


F (j,k) = 


(25) 


where j,k = 1, 


,N b . 


First, we remove the constraint in (21) by introducing a penalty term, 

min ||X|| 2 ,i + ^||^-X —F||| 

xeK N b 2 


(26) 


Here g > 0 is a parameter of the penalty. We will let /x goes to infinity later. Then the 
solution of (26) will converge to the solution of (21). 

Let Y = $0 • X — F, and Y = Y|q = <&q • X — F. Using the Augmented Lagrange 
Multiplier method (ALM) to solve the unconstrained optimization problem (26), we get 
following algorithm. Let Q° = 0, and repeat following two steps until converge 
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( X fc+i Y fc+1 ) = arg min ||X || 2 i + - II'Y " 2 
x,YeM w b xM ’ 2 


2 + |||Y-^-X + F + QV7lll 


• Q fc +! = Q fc + 7 • (Y fc+1 - • X fc+1 + F) 

Here 7 > 0 is a parameter. In above ALM iteration, the main computational load is to solve 
the optimization problem in the first step. Notice that the objective functional depends on 
two terms, X and Y. It is natural to minimize the functional alternately by optimizing one 
of these two terms when the other one is fixed like that in the split Bregman iteration 1113- 


Using this idea, we get following iterative algorithm to solve (26). 
Let Q° = 0 , Y° = 0 and repeat 


X fc+1 = arg min ||X|| 2 ,i + ^\\Y k - $ e • X + F + QV7II2 

xe» N b xM 2 


7, 


Y fc+1 = arg min —IIY||| + — ||Y - $ e • X fc+1 + F 

6 YeR Jv i> xM 2 112 2" 


Q fc /7lll 


Qfc+l = Qfc + 7 . (Y k+1 - . X fc+1 + F) 


fc +1 


r fc+l 


The good news is that the optimization problems in the first and second step can be solved 
explicitly and we can use fast Fourier transform to accelerate the computation. First, let 
us see the first optimization problem. Using the fact that is an orthonormal matrix, we 
have 


X fc+1 = arg min ||X || 2>1 + A||X - 

xeM Jv t xM 2 


7, 


Y k + F + Q fc 


=5 7 ®* e ■ (Y k + F + Q' 


(27) 


where is the conjugate transpose of $ 0 . Here is a shrinkage operator which is defined 
as following. For any v G 


S 7 (v) = 


11 v|| 2—7 

|| v || 2 

0 , 


v, ||v|| 2 > 7, 

||y ||2 < 7 . 


(28) 


If V is a matrix and Vj is its jth row, then S- f (V) is a matrix of the same size as V and the 
jth row is defined to be 5 7 (vj). Notice that is the matrix consisting of Fourier basis 
in 9 coordinate, so the matrix-matrix multiplication can be evaluated by applying discrete 
Fourier transform of each column of the second matrix which can be accelerated by FFT. 


In the second optimization problem, recall that fi > 0 is a penalty parameter in (26) 
which is the larger the better. Now, let fi goes to infinity, then we can solve 


rk+1 


(j, k) = 


($ e ■ X fc+1 — F — Q fc / 7 ) (J, k), (j, k) € n, 
0 , otherwise. 


(29) 


Summarizing above derivation, we get a fast algorithm to solve (21), see Algorithm |4j 
By combining Algorithm [4] and the algorithm in previous section, we can get a complete 
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Algorithm 4 (Fourier extension by group sparsity) 

Input: Q° = 0, Y° = 0. 

Output: Fourier coefficients on over-complete Fourier basis X. 

1: while ||Q a ' — Q fc_1 || 2 > eo do 

2 : X fc+1 = <S 7 -i (T c (Y k + F + Q fc /7)), where T c means apply Fourier transform on 

each column. 

3 , Y k+1 (j, fc) — / (*«-x‘ + 1 -r-Q‘/7)a,*), 

1 0 , otherwise. 

4: Q k+1 = Q k + 7 • {Y k+1 - ■ X fc+1 + F) 

5: end while 

6: X = X fc+1 . 


algorithm to deal with nonperiodic data. Before giving the complete algorithm, there is 
still one issue we need to address. In the derivation of Algorithm [4j we assume that the 
sample points are uniformly distributed in 9 coordinate with given 6. For general signal 
which may not be uniformly sampled in 6 coordinate, we need to do interpolation before 
applying Algorithm |4j 

In order to do the interpolation, first, we need to give the interpolation points which are 
uniformly distributed in 9 coordinate. In this paper, the interpolation points are chosen to 
be 0(0) + (j - 1) A 9, j = 1, • • • , N b , A9 = 2itL e /N b . We set L e = 2 L e and L e = L e(1) 2 ~ e(0) ]. 
This choice is corresponding to the 2-fold over-complete Fourier basis used in ([ 6 ]) . N b = 2 N s 
where N s is the number of sample points of the original signal. The original signals are 
interpolated over the interpolation points by cubic spline interpolation. 

Now, we have an algorithm for nonperiodic signal summarized in Algorithm [5} 


5 Signals with outliers or missing samples 


To deal with the signal with outliers, we need to enlarge the dictionary to include all the 
impulses 

In this case, the optimization problem is formulated in the following way, 


min ||X ||2 l + || Z|| i, subject to: 4 > 0 -X + Z = F. (31) 

xec N t> xM ’ 

ZgRJVsXM 


where and F are same as those in (21) 

Following the similar derivation in Section [4j we obtain Algorithm [ 6 ] to solve above 
optimization problem ( |31| ). By using the same interpolation procedure in Section |4j we can 
integrate Algorithm [ 6 ] with the algorithm in Section [3] to get the method to deal with signals 
with outliers. 

For the signals with missing samples, first, we assign the value of the signal as the 
average of the signal at the locations where the sample was lost and then treat the missing 
sample as the outliers. 
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Algorithm 5 (Algorithm for nonperiodic signal) 


Input: Initial guess of the phase function 9 k = 9q, ?] = 0. 

Output: Phase functions and the corresponding envelopes: 9k, al, j = !,■■■ ,M. 


2 

3 

4: 

5: 

6 : 

7: 


while 77 < A do 

while || 9 k +1 — 0£|| Z 2 > e 0 do 

Interpolate 77 . _1 to a uniform mesh in the ^((-coordinate to get rgn. 

Using Algorithm|4]to get the Fourier coefficients of rgn on the over-complete Fourier 
basis. 

Apply a cutoff function to the Fourier Transform of rgn to compute a and b on the 
mesh of the ^-coordinate, denoted by ag « and bgn. 

Interpolate ag and bg » back to the uniform mesh of t. 

Calculate averaged update of frequency: 


8 : 


Aujj = 


T +1 (e + 1 )'-r +1 (°r +1 )' 

•f + 1 ) 2 +(e +1 ) 2 


Au; = 


xv/V 
2 ^ 7=1 1 k 


(30) 


where r £ n+1 = (< 4 ,Tl+1 ) + and (•) / denote the derivative of (•) with 

respect to t. 

Update 6% 

A9'= P v{ g. v) (Aco), A9(t)= f A9\s)ds, 9™ +1 = 9 n k - (3A9, 

Jo 

where /? G [ 0 , 1 ] is chosen to make sure that 9 k +1 is monotonically increasing: 

/3 = max j cc G [0,1] : — {9 k — aA9) > 0 


and Pv p (6-,g) is the projection operator to the space V p (9] rj). 

9: end while 

10 : 77 = 77 + Af] 

11: end while 
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Algorithm 6 (Algorithm for signals with outliers) 

Input: Q° = 0, Y° = 0. 

Output: Fourier coefficients on over-complete Fourier basis X, estimate of outliers Z. 

1: while || Q a ' — Q fc_1 || 2 > eo do 

2: X fc+1 = 5,,-i (j r c ( Y k — Z k + F + QVy))> where T c means apply Fourier transform 

on each column. 

3, Y‘+‘(i,fc) = { (^■X l + 1 + Z t -F-QV 7 )a* ; ), (j,k)eSi, 

' 0, otherwise. 

4: Z k+1 = Ty-i (Y k+1 + F + Q fc /7 — • X), where is a shrinkage operator. 


For any i£l, 


'Fy(x) 


M~7 

M 

0 , 


x\ > 7 , 
x\ < 7 . 


5 

6 
7 


For any matrix A = ( ajk ), 7iy(A) = (7iy(ajfe)). 
Qfc+i = Qfc + 7 . (z*+i _ . X fc+1 + F) 

end while 


X = X 


fc+i 


z = z 


k +1 


o, is defined in (24). 


6 Numerical results 

In this section, we use two numerical examples, one is synthetic and one is real data, to 
demonstrate the performance of the method proposed in this paper. The first example is a 
simple synthetic signal which is used to test the robustness of the method to the perturbation 
of white noise. 

Example 1: Synthetic signal with white noise 

In this example, the signals are generated as following: 

f(t) = cos(40vr(f + l) 2 ), f(t) = f(t) + 5X J (t), j = 1, • • • , 10, t G [0,1]. (32) 

where X 3 (f), j = 1, • • • ,10 are independent white noise with standard deviation a 2 = 1. 
The number of sample points is 512 and the sample points are uniformly distributed over 

[o,i]. 

The original signals are plotted in the left figure of Fig. [l} As we can see that the noise 
is so large such that we can not see any pattern of the original clean signal. The recovered 
instantaneous frequency is shown in the right figure of Fig. [l] If we recover the frequency 
from each signal separately, since the noise is too large, the frequencies are totally wrong. 
However, if we use the structure that 10 signals have same instantaneous frequency, the 
frequency recovered is much better. 

Example 2: Cable tension estimat^] 

1 The authors would like to thank Prof. Yuequan Bao of the School of Civil Engineering of Harbin 

Institute of Technology for the experimental data and for permission to use it in this paper. 
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Figure 1: Left: 10 measurements generated by f(t) + 5 X{t). Bold black: signal with¬ 
out noise; thin lines: 10 measurements with larger noise; Right: instantaneous frequency 
recoverd from signal with 10 different measurements with larger noise. Bold blue: exact 
frequency; bold black: recovered frequency from 10 signals; thin lines: frequency obtained 
from 10 signals separately. 

This example is a real problem about the tension estimate of the cables in bridge. Before 
demonstrating the numerical results, we give some backgrounds first. 

For large span bridges, such as cable-stayed bridge and suspension bridge, the cables 
are a crucial element for overall structural safety. Due to the moving vehicles and other 
environmental effects, the cable tension forces vary over time. This variation in cable tension 
forces may cause fatigue damage. Therefore, estimation of the time-varying cable tension 
forces is important for the maintenance and safety assessment of cable-based bridges. 

One often used method to estimate the cable tension force is based on the instantaneous 
frequency estimate of the cable vibration signal. According to the flat taut string theory 
that neglects both sag-extensibility and bending stiffness, cable tension force, F, can be 
calculated by 

m = imL2 {^) 2 (33 > 
where uj n (t) is the time-varying nth natural frequency in radius/s and m, L are mass 
density and length of cable. 

Also from the flat taut string theory, an important and useful feature of the vibriations 
of the cable is that the natural frequencies of the higher modes are integer multiples of the 
fundamental frequency, that is oj n (t) = nui(t). This feature means that we can combine the 
information of different modes together to recover the instantaneous frequency. Then the 
method developed in this paper for multiple signals applies after some minor modifications. 
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Figure 2: Left: vibration signal of the cable; Right: zoom in of the signal to demonstrate 
the missing samples (flat segments). 

The original experimental signal is given in Fig. [2] Obviously, the signal has some 
outliers. In the original signal, about 10% of the samples were lost which is demonstrated 
very clearly in the zoom in picture in Fig. [2j The tension force estimated by our method 
is given in Fig. [3] and [4j The tension forces obtained from different modes are shown in 
Fig. [3j If only one mode is used, the estimation of the force is not very accurate. There 
are many oscillations although the overall picture match with the measured force. If we use 
1-5 modes together to estimate the tension force, the result is much better, Fig. [4j 

7 Conclusion 

In this paper, we consider multiple signals and these signals share the same instantaneous 
frequencies. This kind of signals emerge in many scientific and engineering problems. By 
exploiting the structure of the common frequencies, we develop some algorithms to find the 
sparsest time-frequency decomposition. These algorithms can be accelerated by FFT, so 
they are very efficient. The algorithms are also very robust to noise, since the information 
of multiple signals are used simultaneously. 
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Figure 3: 


Cable force estimated by single natural mode 1-5 from up to bottom. 
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Figure 4: Cable force estimated by using 1-5 natural modes together. 
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